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ABSTRACT 

We study the stability of mean-motion resonances (MMR) between two planets during their migration in a protoplanetary disk. We 
use an analytical model of resonances, and describe the effect of the disk by a migration timescale (T^j) and an eccentricity damping 
timescale (T^j) for each planet (I = 1,2 respectively for the inner and outer planet). We show that the resonant configuration is stable 
if Tg i/rg 2 > (^ 1 /^ 2 )^- This general result can be used to put constraints on specific models of disk-planet interactions. For instance, 
using classical prescriptions for type I migration, we show that when the angular momentum deficit (AMD) of the inner orbit is larger 
than the outer’s orbit AMD, resonant systems must have a locally inverted disk density profile to stay locked in resonance during the 
migration. This inversion is very untypical of type I migration and our criterion can thus provide an evidence against classical type 
I migration. That is indeed the case for the Jupiter-mass resonant systems HD 60532b, c (3:1 MMR), GJ 876b, c (2:1 MMR), and 
HD 45364b, c (3:2 MMR). This result may be an evidence for type II migration (gap opening planets), which is compatible with the 
large masses of these planets. 
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1. Introduction 


In |Delisle et al.| ( [20T^ , we showed that tidal dissipation raised 
by the star on two resonant planets can produce three kinds of 
distinct evolutions depending on the relative strength of the dis¬ 
sipation in both planets. The three different outcomes of this tidal 
process are systems that stay in resonance, systems that leave 
the resonance with an increasing period ratio {Pout IP in) and sys¬ 
tems that leave the resonance with a decreasing period ratio. For 
known near resonant systems, the comparison of the period ratio 
of the planets with respect to the nominal resonant value helps to 
put constraints on the tidal dissipation undergone by each planet 
and thus on the nature of the planets (see |Delisle et al.|2014| ). In 
this article, we generalize our reasoning to other forms of dissi¬ 
pation, in particular to disk-planet interactions. 


Disk -planet interactions can induc e migration of the plan¬ 
ets (e.g. [Goldreich & Tremaine||1979] ). In the case of conver¬ 
gent migration (i.e. decreasing period ratio), the planets can be 
locked in resonance (e.g. [Weidenschilling & Davis||1985| ). Two 
planets that are locked in resonance have their eccentricities ex¬ 
cited on the migration timescale (e.g. [Weidenschilling & Davis | 
|1985| ). However, disk-planet interactions also induce exponential 
eccentricity damping. Depending on the respective timescales of 
the migration and eccentricity damping, the system can reach 
a stationary state in which eccentricities stay constant ( [Lee &| 
Peale|[2005] ). The semi-major axes continue to evolve but the 


semi-major axis ratio (or period ratio) stays locked at the res¬ 
onant value. Recently, [Goldreich & Schlichting ( [2014[ ) showed 
that this equilibrium is unstable in the case of the circular re¬ 
stricted three body problem where the inner planet has negligible 
mass. This means that after the resonance locking, the eccentric¬ 


ity of the inner planet reach an equilibrium value but then under¬ 
goes larger and larger oscillations around this equilibrium value 
until the system reaches the resonance separatrix and leaves the 
resonance. Then, the period ratio is no more locked at the reso¬ 
nant value and the convergent migration continues (decreasing 
period ratio) until the system reaches another resonance. The 
timescale of the resonance escape is given by the eccentricity 
damping timescale and is thus short compared to the migration 


[dreich & Schlichting[ ( [2Q14[ ) concluded that when the disk dis¬ 
appears and the migration stops, only a few systems should be 
observed in resonance. However, this conclusion is mainly based 
on a particular case in which the mass of the inner planet is much 
smaller than the mass of the outer planet whose eccentricity is 
negligible and the migration and damp ing forces are only un der- 
gone by the inner planet. As shown in [Delisle et al.[ ( [MT^ , the 
evolution of a resonant system under dissipation highly depends 
on which planet is affected by the dissipation. In this paper we 
study a more general case in which both planets have masses, 
eccentricities, and undergo dissipative forces. 

In Sect. [ 2 ] we introduce the notations and the model of the 
resonant motion in the conservative case that we developed in 
[Delisle et al.[ ( 2014} . In Sect.j^we study the dissipative evolution 
of resonant planets in a very general framework (Sect. [3.1[), and 
we apply this modeling to disk-planet interactions (Sect. [3.2| ). In 
Sect. we show how our model can be used to put constraints 
on disks properties for observed resonant systems. In Sect. 
we apply these analytical constraints to selected examples, and 
compare them with numerical simulations. We specifically study 
HD 60532b , c ( 3:1 resonance. Sect. [5.1[ ) GJ 876b, c (2:1 reso¬ 
nance, Sect. |5J, and HD 45364b, c (3:2 resonance, Sect. [5.3[ ). 


timescale (see [Goldreich & Schlichtin^[2014[). Therefore, [Gol- 
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2. Resonant motion in the conservative case 

In the following, we refer to the star as body 0, to the inner planet 
as body 1, and to the outer planet as body 2. We note the 
masses of the three bodies, and introduce pt = Q(mo + mt) and 
Pi = moMilimo + nii), where @ is the gravitational constant. We 
only consider the planar case in this study. 

In [Delisle et al.| ( |2014| ) we constructed a simplified, and in- 
tegrable model of the resonant motion in the conservative and 
planar case. The main simplification of this model is to assume 
that the eccentricity ratio (ei /^ 2 ) stays close to the forced eccen¬ 
tricities ratio (e I ^eii 1^2,611)- These forced eccentricities correspond 
to the eccentricities at the elliptical fixed point at the resonance 
libration center. With this assumption, and assuming moderate 
eccentric ities p] the H amiltonian of the system can be simpli¬ 
fied (see [Delisle et al.||2014| ) to the following simple pendulum 
Hamiltonian 


^ = -6^ cos(q6), 


( 1 ) 


where q is the degree of the resonance (q = k 2 - h for a ^ 2 -^i 
resonance), e is the action coordinate and provides a measure of 
the distance to the exact commensurability. 6 is the unique res¬ 
onant angle in this simplified model. It is a combination of both 
usual resonant angles (ct/ = ^42 - see Appendix [Bj 

Eqs. ( |B.4| ) and ( |B.6| ) and [Delisle et al.||2014| ). R is sl constant 
that depends on the masses of the bodies and on the considered 
resonance (see jDelisle et al.|2014| ). d is a constant of motion (pa¬ 
rameter of the model). We have 


^ = Ai - Ai^o + A2 - A2,o, 

^ = ^1,0 - Gi A 2 ,o - G2, 


( 2 ) 

( 3 ) 


where A/ is the renormalized circular angular momentum of 
planet /, Gt its renormalized a ngular momentum (see Ap¬ 
pendix]^ and [DeIIsl££t^^|20^. The subscript 0 denotes the 
values at the exact commensurability. The quantities A/ only de¬ 
pend on the semi-major axis ratio a = a\la 2 


Ai(a) = 


A2((r) = 


1 


(^2 Ai) + iPiiPx) Vi^wo^) 

1 

{k2lki) + (m2/mi)/V^’ 

1 


1 -r {k 2 lki){PxlP 2 ) ^Jpla|p 2 
1 

1 -I- (^ 2/^1 )(^i 7 ^ 2 ) 

At the exact commensurability, we have 


(4) 


(5) 




ki \ 


2/3 


( 6 ) 


The quantities Gi depend on a and on the planet eccentricities 
Gi(a, e,) = A, (O') -e]. 


(7) 


^ The pendulum approximation of resonances is obtained using an an¬ 
alytical expansion in power series of eccentricities and is thus not valid 
at high eccentricities. Moreover, when eccentricities are vanishing, the 
phase space bifurcate and a better approxi mation is given by the sec¬ 
ond fundamental m odel of resononaces (see jHenrard & Lemaitre| 19831 
Delisle et al.|2012| l. 
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Fig. 1. Phase space of a resonance of order q in the simplified pendulum¬ 
like approximation (Hamiltonian 6 is the unique resonant angle 
and 6 its conjugated action. The separatrix is highlighted in red. The 
amplitude A (defined with O^ax see Eq. |T^ ) is 0 at the center of the 
resonance (elliptical fixed point) and 1 at the separatrix. 


Let us note p the renormalized angular momentum deficit 
(AMD) of planet i ( |Laskar|1997|[2000l i 

/,• = A,- - Gi = ^Ai^f oc el (8) 


with 




(9) 


The simplifying assumption introduced in Delisle et al. (20141 
implies (see also Appendix |B]) 

h h,eii . 2 , 

- = 7 — = tan^ 0 , 

h h,ell 


( 10 ) 


where 0 is a constant angle and l^eii are values of the renormal¬ 
ized AMD at the center of the resonance (elliptical fixed point, 
seejDelisle et al.||2014j). We also note D the renormalized total 
AMD 


!D = li+l2 = 6 + e. ( 11 ) 

The parameter 6 corresponds to the renormalized total AMD at 
the exact commensurability (S = Dq). Thus, for a resonant sys¬ 
tem, S provides a measure of the planet eccentricities (S oc 
see Eq. ^). Eigure shows the phase space corresponding to 
Hamiltonian m. The width of the resonant area is proportional 
to oc for a resonance of order q (see Eig.[^. Eor a reso¬ 
nant system, in the regime of moderate eccentricities, a measure 
(between 0 and 1) of the relative amplitude of libration (ampli¬ 
tude of libration versus resonance width) is given by (seejDelislej 
jet al.|2014j ) 

A = sin2(^), (12) 

where Omax is the maximum value reached by the resonant angle 
6 during a libration period (see Eig.[^. 

Note that our simplifying assumption (eccentricity ratio 
close to the forced eccentricities ratio) is well verified when the 
amplitude of libration is small (A 1) and the system stays 
close to the elliptical fixed point. Eor high amplitude of libration 
(A ~ 1), the eccentricity ratio undergoes oscillations around the 
forced value and our model only provide a first approximation 
of the motion (see jDelisle et al.|2014j ). 
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3. Resonant motion in the dissipative case 

In this section we describe the evolution of a resonant system 
undergoing dissipation. The main parameters that have to be 
tracked during this evolution are the parameter S which describes 
the evolution of the phase space (and of the eccentricities for res¬ 
onant systems) and the relative amplitude A which describes the 
spiraling of the trajectory with respect to the separatrix of the 
resonance. 

3.1. General case 


where Gt is the angular momentum of planet i. From these sim¬ 
ple decay laws we can deduce the evolution of t he pa rameters of 
interest for resonant systems (S and A, see Sect. |3.l| ). We have 



(19) 


at 

d 


+ 2- 




_ t 



( 20 ) 


Let us consider a dissipative force acting on the semi-major axes 
and the eccentricities of both planets. We first consider a very 
general case and do not assume a particular form for this dissi¬ 
pation, except that it acts on a long timescale. The evolution of 
the system can be described by the three following timescales 
(which may depend on the eccentricities and semi-major axes of 
the planets): (^ 2 /^ 2 )^/, (a/a)d- Note that, for sufficiently 

small eccentricities, we have and 


(13) 


The evolution of the parameter S that drives the evolution of 
the phase space (and of the eccentricities for resonant systems) 
is given by (see Appendix 


5\d = 2(cos^^^ 


+ sin^ (j) — 


6 


D 


A 2 - sin^ (p a 
2 a 
q A 1 A 2 a 
ki 2 a 


D 


(14) 


: A >= - (< €€\d > < €^S\d >] , 

2RS^/^ \ 4S I 


with 


c\d 


q A 1 A 2 a 
k\ 2 a 


(15) 


(16) 


3.2. Disk-planet interactions 

Let us now apply Eqs. ([Tg, ([g to the specific case of disk- 
planet interactions. Because of these interactions the planets un¬ 
dergo a torque that induces a modification in their orbital el- 
ements and subsequent migration in the disk (e.g. [Goldreich 


|paloizou|2007t|Goldreich & Schlichting|2bl4|): 


Gi 

ei_ 


1 

Tm,i ’ 

1 


(17) 

(18) 


The evolution of the semi-major axis ratio is thus governed by 


1 

to 

h 

^ Ai sin^ 0 

A 2 COS^ (f)\ 

d Tm A 1 A 2 

\ Tea 

Te,l ] 


with 

1 

Tm 


1 

Tm,2 


1 

Tm,l 




( 21 ) 


( 22 ) 


From Eq. ( p^ we obtain 
q A 1 A 2 




h 

, A2-sinV ^ 

Tm 

^ I cos^ (b sin^ 0 \ ^ 


Te,\ Te,2 
q ( Ai sin^ 0 Ai cos^ 0 




For a resonant system, the evolution of the relative amplitude of 
libration reads (see |Delisle et al.|2014[ Appendix A) 


-r2- , 

kl \ Te,2 

q A 1 A 2 

kl Tm 

-D 


D 


(23) 


^. s. ^2 cos^ 0 sin^ 0 
2(Ai+A2)It^^:^+ ^ 


kl Te, 


G,2 


A 2 - sin^ 0 


where we neglect second order terms in D oc e^). Let us 
note 


|& Tremam^|1979| |1980| ). In particular, the angular momentum 
of each planet evolves on an exponential timescale Tm,i due 
to this migration, while eccentricities evolve on an exponential 
timescale (e.g.|Papaloizou & Larwood|2000[ Terquem & Pa-| 


1 q A 1 A 2 

Tm kl Tm 

1 ^.. . . / ^2 cos^ 0 sin^ 0 \ Ao- sin^ 0 

— = 2(Ai+A2) + ^ 

J^E \M ^6,2) J-m 

We thus have 

AI ^ ^ 

0\d = — -—, 

f M dE 

<^\d> = 


(24) 

(25) 

(26) 
(27) 


Note that the damping timescale is often much shorter than the 


migration timescale (Tei Tmi, c.g. Goldreich & Tremaine 
[T9§01 ), thus 


1 ^. s / ^2 cos^ 0 sin^ 0 
— 2(Ai -r A 2 ) -—-— + —— 

dE \M le,2 


(28) 
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The timescales Te,Tm can be expressed using more usual nota¬ 
tions 


1 

Tm 


1 

Te 


q ( k2 m2 1 IkA^ nil /—I * 

ki [ ki nil 1^2/ m 2 

1 1 

Tm,2 Tjn,\ , 

2fl + !!Ilv55)fi4=lV5;r'x 

\ m 2 / \ ^1 m 2 / 


Te,2 kimVi), 


/ell 


^e,\ 


m 2 \^2 


(29) 


(30) 


x-l 




ell 


Let us now compute the evolution of this amplitude of libra- 
tion. According to Eq. ([Tg, we need to compute < ee\d > and 
< e^d\d >. We have 

^ -I .. g /A2COSV Aisin^ 0 \ 2^ 

< ee\a > = 2— —---- < e >, (32) 

h \ Te,l Te,2 I 

<£^d\d> = >=<5\d><e^ >, (33) 

where < > can be computed using elliptic integrals (see 

IDelisle et al.lMT^ 


< e > 


IRS^'^A. 


(34) 


Depending on the values of Tm and Te, different evolution 
scenarios for S are possible. Note that all these equations remain 
valid for Tm and Te negative. In most cases the disk induces a 
dampin g of eccentricities {Te^ > 0, thus > 0), but some stud¬ 
ies (e.g. Goldreich & Sari|2003| ) suggest that an excitation of the 
eccentricities by the disk is possible (Tg / < 0, thus Te < 0). 
Tm is positive if the period ratio between the planets {Pi!Pi) 
decreases (convergent migration). But if the planets undergo di¬ 
vergent migration (^ 2/^1 increases), is negative. This does 
not depends on the absolute direction (inward or outward) of the 
migration of the planets in the disk but only on the evolution of 
their period ratio. 

In the case of divergent migration, the planets cannot get 
trapped in resonance (e.g. [Henrard & Lemaitre|1983| ). The sys¬ 
tem always ends-up with a period ratio higher than the resonant 
value and this does not depend on the damping/excitation of ec¬ 
centricities. 

The case of convergent migration is more interesting. If the 
initial period ratio is higher than the resonant value, the plan¬ 
ets can be locked in resonance. This induces an excitation of the 
eccentricities of the planets {d\M = ^ITm > 0). If Te <0 (excita¬ 
tion of eccentricities by the disk) or » Tm (inefficient damp¬ 
ing), S (as well as the eccentricities) does not stop increasing. 
When eccentricities reach too high values, the system becomes 
unstable and the resonant configuration is broken. 

The most common scenario is the case of efficient damping 
of eccentricities (0 < < T^). In this case, S reaches an equi¬ 

librium value (< S\d >= 0, see Eq. (27)) 


^eq = 


h. 

Tm 

q A 1 A 2 / ^ _^ 

2 ki Ai -r A2 \ Tyn,2 Tm,i 

\2 \ 

X 

/ell 


2 ^\“ 


k 2 cos^ 0 ^ sin 0 


kl Te,l 


f e,2 


2 1 m 2 \e2 


7 7 7 , m 2 I 

^1+^2+ ^2- yOlQ kl - - - 

m 2 mi 


1 

Tm,2 


1 


* m,l 


Te,2 kim2\e2]^ii Te 


Note that the first term (< ee\d >) does not depend on the mi¬ 
gration timescale but only on the damping timescale. The second 
term (< e^6\d >) vanishes when the system reaches the equilib¬ 
rium d = deq, since d\d = 0. This is not surprising because the 
first term describes the evolution of the absolute amplitude of li- 
bration 6^, while the second one describes the evolution of the 
resonance width which does not evolve if the phase space does 
not evolve (constant d). Einally, we obtain (see Eq. 


(A 2 cos^ 0 Ai sin^ 0 


kl 


^ e,l 


f e,2 


The amplitude of libration increases if 


A 2 cos^ 0 ^ Ai sin^ 0 




^ e,2 


which is equivalent to 

\2 


Li< 

E,2 W 2 


ell 


Using ^ eu this gives 


A/ 


TeA ^ 

Te,2 \^2 


(35) 


(36) 


(37) 


(38) 


ell 


(31) 


However, as shown by [Goldreich & Schlichting| ( |2014| ) for the 
restricted three body problem, this equilibrium can be unstable. 
Even if the parameter 6 reaches the equilibrium 6eq and the phase 
space of the system stops evolving, the amplitude of libration can 
increase until the system crosses the separatrix and escapes from 


resonance. 


where the eccentricity ratio is evaluated at the elliptical fixed 
point {ell subscript) at the center of the resonance. In the circu¬ 
lar restricted case studied by [Goldreich & Schlichtingj ( [2Q14[ ), 
^2 = 0 and Te ;2 - thus the amplitude always increases 
(Eq. ( [38] )) and the equilibrium is unstable. However, in the op¬ 
posite restricted case {ei - 0 and Tqa - + 00 ), that was not ad¬ 
dressed in [Goldreich & Schlichting[p014[ ), the amplitude always 
decreases (Eq. ([38[)) leading to a stable equilibrium. 

Note that this result is based on our approximation that the 
eccentricity ratio remains close to the forced value (value at the 
elliptical fixed point). This is well verified for a small amplitude 
of libration but when the amplitude increases, the eccentricity ra¬ 
tio oscillates and may differ significantly from the forced value. 
Our model thus only provides a first approximation of the mean 
value of this ratio in the case of high amplitude of libration. 

To sum up, the evolution of a resonant pair of planets un¬ 
dergoing disk-planet interactions depends mainly on two param¬ 
eters: TeITm (damping vs migration timescale) and 1 / 7^2 
(damping in inner planet vs outer planet). The ratio TeITm gov¬ 
erns the equilibrium eccentricities of the planets (see Eqs. S’ 
([3T])). The ratio Tg 1 / 7^2 governs the stability of this equilibrium 
(see Eqs. 
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4. Constraints on disk properties 


In this section, we show how the classification of the outcome of 
disk-planet interactions can be used to put constraints on the dis¬ 
sipative forces undergone by the planets and thus on some disk 
properties. More precisely, if a system is currently observed to 
harbour two planets locked in a MMR, it is probable that this 
resonant configuration was stable (or unstable but with a very 
long timescale) when the disk was present. We could imagine 
that the configuration was highly unstable but the protoplane¬ 
tary disk disappeared before the system had time to escape from 
resonance, however this would require a fine tuning of the disk 
disappearing timing. Thus the amplitude of libration was prob¬ 
ably either decreasing, or increasing on a very long timescale. 
This induces that 




(39) 


Moreover, a small amplitude of libration is probably the sign 
of a damping of the amplitude on a short timescale 


^»(-f 

Te,2 \ ^2 / ell 


(40) 


model (see Eq. ( [ST] )) which themselves depend on the proper¬ 
ties of the disk and the planets. There exists a wide diversity of 
disk models in the literature, which would result in significantly 
different migration and damping timescales for each planet. Our 
analytical model is very general and can handle these different 
models as long as expressions for i, 7^ 2 , and are 
available. 

We consider here the case of type I migration to illustrate 
the possibility of constraining the disk properties for observed 
systems. Following the prescriptions of |Kley & Nelson] ( |2012| ), 
we have 


Tm,\ _ m 

Tm,2 OTi V fli \//(a2)/a2/ 

Zk ~ iHjad? 

Tm,i \ / 


(44) 

(45) 


where H{a) is the local disk scale height and lL{d) oc its 
local surface density. The standard MMSN (Minimum Mass So¬ 
lar Nebula) model assumes E oc = 3/2). HI a is called 

the disk aspect ratio and is o ften assumed to be rou ghly constant 
and of the order of 0.05 (e.g. |Kley & Nelson|2012| ). Using these 
assumptions we obtain 


On the opposite, a large amplitude of libration could be the sign 
of a long timescale of amplitude damping or a long timescale 
of amplitude excitation. Indeed, if the amplitude was increasing 
fast, the system should not be observed in resonance. If it was 
decreasing fast, the observed amplitude should be very small. 
However, another mechanism may be responsible for the excita¬ 
tion of the amplitude of libration, possibly after the disk disap¬ 
pearing (e.g. presence of a third planet in the system). Thus we 
cannot exclude the case of a fast damping of the amplitude of 
libration, even in the case of an observed large amplitude. 




(41) 


In addition to the constraints obtained from the observed am¬ 
plitude of libration, the observed values of both eccentricities 
is also an important information. If the planets did not undergo 
other source of dissipation since the disk has disappeared, the 
present eccentricities should still correspond to the equilibrium 
ones. For close-in planets, the tides raised by the star on the plan¬ 
ets induce a significant dissipative evolution of the system after 
the disk disappearing (e.g. [Delisle & Laskar|[2014| ). Therefore, 
this reasoning only applies for systems farther from the star for 
which tidal interactions have a negligible effect on the orbits over 
the age of the system. Let us recall that the equilibrium eccen¬ 
tricities are given by (Eq. ( [ST] )) 


Tm,2 

Tm,l _ ^m,2 IH \~^ 

Te,l Te,2 

For the sake of brevity, we note in the following 


(46) 

(47) 


^2 ^2-1/2 ^ ^ 

^1 Tm,2 Te,2 


^ Tm,2 

V <2 / Te l 7^ 2 


(48) 

(49) 


If the system is observed with a small amplitude of libration in 
the resonance, we have (Eq. ( [4Q| )) 


T » 




\^2/ell 

and if the amplitude is large we have (Eq. <ED) 

\2 


r>|^ 

^2 


(50) 


(51) 


ell 


The lower limit we obtain for r corresponds to an upper limit 
for the density profile exponent (see Eq. ([48])). In particular, 
if 


S = Seq = 


II 

Tm 


with 


1 


-(Ai^?+A2^2') 


k2 ^ m2 1 

h Ml ^/a 


^1 


k2mi j- 

1 -r -- Sa 

ki m2 


_2. 

2 ■ 


(42) 

^1 


(52) 


(43) 


the density profile of the disk must be inverted < 0, i.e. the 
surface density increases with the distance to the star) for the 
system to be stable in resonance. The condition of Eq. is 
roughly equivalent to 

h > /2, (53) 


S can be computed from the known (observed) orbital elements 
of the planets. Thus, the ratio TeITm of the damping and mi¬ 
gration timescales is constrained by the observations. This ratio 
depends on the four timescales (T^j, 7^,2, Tm,\, and Tm, 2 ) of the 


where h , h are the angular momentum deficits (AMD) of both 
planets. 

Let us recall that in order to be captured in resonance the 


planets must undergo convergent migration (e.g. Henrard & 
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[Lemaitre 


(see Eq. (48|)) 

T> 1. 


1983]). This put another constraint on the parameter r 

iwr 

(54) 


Again, this corresponds to an upper limit for /3-z, and if 


< 1 , 

mi 


(55) 


the density profile of the disk must be inverted < 0) for 
the planets to undergo convergent migration. The condition of 
Eq. ( [55] ) is roughly equivalent to 


Ai > A 2 , 


(56) 


where Ai, A 2 are the circular angular momenta of both planets. 

The constraint provided by the observation of the equilib¬ 
rium eccentricities reads (see Eq. 


Cl T-l 


K T + C 2 


(57) 


with 

Cl 




^ ^2 m\ , _ k\m2 1 

1-1--I-— 

q m 2 q mi 


C 2 


k2 mi I ei 


ki m 2 \e2 


v^, 


ell 


and deq is given by Eq. (43). We thus have 

Cl T-l 


K = 


S T + C2' 


(58) 


(59) 


(60) 


where Ci, C 2 and d can all be derived from the observations. 
Note that K is an increasing function of r (Eq. (|60|)). Thus, our 
analytical criterion for stability provides a lower bound for both 
T and K. 


5. Application to observed resonant systems 

In the following we apply our analytical criteria to systems that 
are observed in resonance. We also performed N-body simula¬ 
tions with the additional migration and damping forces exerted 
by the disk on the planets. We used the ODEX integrator (e.g 
Hairer etar]|2010 ) and the dissipative timescales (angular 
momentum evolution), (eccentricity evolution) are fixed for 
each planet and each simulation. 

Many multi-planetary systems are observed close to differ¬ 
ent MMR (period ratio close to a resonant value). However only 
a few of them have a precise enough determination of the plan¬ 
ets orbital parameters in order to distinguish between resonant 
motion and near-resonant motion. To our knowledge, all known 
resonant planet pairs are giant planets (better precision of orbital 
parameters). We thus selected three of these resonant giant planet 
pairs to illustrate our model. 

Note that giant planets are believed to undergo type II mi¬ 
gration. Our analytical model is very general and can take into 
account any prescription for the evolution of the angular momen¬ 
tum and the eccentricity of each planet. We did not find a sim¬ 
ple analytical prescription for type II migration in the literature. 
Indeed type II migration is a more complex (non-linear) mech¬ 
anism than type I migration, and it is not yet well understood. 
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Tabl e 1, Orbital parameters o f HD 60532b,c used in this study (taken 
from |Laskar & Correia|2009| ). The stellar mass is 1.44Mo. 


Parameter 

[unit] 

b 

c 

m 

[Mj] 

3.1548 

7.4634 

P 

[day] 

201.83 

607.06 

a 

[AU] 

0.7606 

1.5854 

e 


0.278 

0.038 


In particular, the timescale of type II migration is still discussed 
(e.g. [Duffell et al.|2014 |Durmann & Kley|2015| ). However, the 
effect of type II migration is expected to be similar to type I mi¬ 
gration (i.e. inward migration and damping of eccentricities, e.g. 
[Bitsch & Kley||2010| ). The main difference is that the disk pro¬ 
file is affected by the presence of the planets (gap around the 
planets) and the timescales of migration and damping are thus 
slowed down. We thus chose to apply type I migration prescrip¬ 
tions for our study of these giant planets resonant systems as a 
first approximation. 


5.1. HD 60532b, c: 3:1 MMR, large amplitude of libration 

The star HD 60532 hosts t wo planets (see|Desort et al.|2008| ) that 
exhibit a 3:1 period ratio. [Laskar & Correia] ( |2009| ) performed a 
dynamical study of the system and confirmed the 3:1 MMR be¬ 
tween the planets with a large amplitude of libration (~ 40°). 
We reproduced the orbital elements of the planets (taken from 
[Laskar & Correia||2009| ) in Table Eor this system, the forced 
eccentricity ratio (ratio of eccentricity at the center of the reso¬ 
nance) is 



( 61 ) 


Since the system is observed with a large amplitude of libration, 
the stability constraint gives (Eq. ( [5T] )) 


T > 



(62) 


Note that this value is much larger than one, thus the condition 
of convergent migration is fulfilled (see Eq. g). This value 
of T corresponds to a surface density profile exponent of (see 
Eq. 


< -1.3. 


(63) 


Let us recall that for the MMSN model Py = 3/2. The negative 
value we obtain corresponds to an inverted density profile. This 
inverted density profile is very untypical of type I migration. This 
result is thus an evidence that the planets did not undergo a clas¬ 
sical type I migration. This is not surprising since giant planets 
are expected to open a gap and undergo type II migration. Our 
results also constrain a type II migration scenario. Indeed, inde¬ 
pendently of the migration prescriptions, for the resonance to be 
stable the damping of the outer eccentricity must be much more 
efficient than the inner eccentricity damping {Te,ilTe ,2 ^ 9). One 
would need prescriptions for type II migration to relate this result 
with some disk and/or planets properties. 

Erom the observed orbital elements we obtain (see Eqs. ( |4^ , 
( [58] ), and ( [59] )) 


6 

= 6.5 X 10-^ 

(64) 

Cl 

= 0.44, 

(65) 

C 2 

= 7.9. 

(66) 
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Fig. 2. Semi-major axes, period ratio, eccentricities, eccentricity ratio, and angles evolution for simulations of HD 60532b, c with different 
dissipation timescale ratios r = The ratio K - is set according to Eq. ( [60| ) to reproduce the observed equilibrium 

eccentricities. We used r = -i-oo, 10, 8, 4 with K = 70, 34, 30, 17 respectively for the four shown sirnulations (four columns). The amplitude of 
libration decreases for the first two simulations (r = -hoo, 10) and increases for the last two (r = 8, 4). The value given by our analytical criterion 
for the transition between decreasing and increasing amplitude is r ~ 9. 


Using T > 9, the current eccentricities should be reproduced with 
(see Eq. ( [60| ) 

K > 30, (67) 

which corresponds to an aspect ratio of (see Eq. ( [49| )) 

- < 0.18. (68) 

a 

We performed numerical simulations with different values 
of T. Eor each simulation, the value of K is computed using 
Eq. ( [60| ), in order to reproduce the equilibrium eccentricities. 
We fixed 7^,2 = 5 x 10^ yr for all simulations, and integrated 
the system for 10^ yr. We thus have 7tn,\ = 5 x lO^r yr, 


7^ 2 = 5 X 10^/7^ yr, i = 5 x lO^r/T^ yr. The semi-major axes 
are initially 10 and 22 AU, such that the system is initially out¬ 
side of the 3:1 resonance with a period ratio of about 3.3. Both 
eccentricities are initially set to 0.001 with anti-aligned perias- 
trons and coplanar orbits. The planets are initially at periastrons 
(zero anomalies). The evolution of the semi-major axes, the pe¬ 
riod ratio, the eccentricities, the eccentricity ratio, and the angles 
are shown in Eig.[2 These simulations confirm our analytical re¬ 
sults: for T < 9, the amplitude of libration decreases, and for 
T > 9 the amplitude increases (see Eig.|^. 
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Fig. 3. Same as Fig. [^but for GJ 876b, c. We used r = +oo, 20, 10, 5 with K = 130, 58, 36, 19 respectively for the four shown simulations (four 
columns). Note that the last simulation (r = 5, = 19) ended before 10^ yr (around 8 x 10^ yr) because of orbital instability when the system 

escaped from resonance. The amplitude of libration decreases for the first two simulations (r = +oo, 20) and increases for the last two (r = 10, 5). 
The value given by our analytical criterion for the transition between decreasing and increasing amplitude is r ~ 42. 


Table 2. Orbital parameters of GJ 876b,c used in this study (taken from 
[Correia et al.|2010| ). The stellar mass is 0.334Mo. 


Parameter 

[unit] 

c 

b 

m 

[Mj] 

0.86 

2.64 

P 

[day] 

30.259 

61.065 

a 

[AU] 

0.132 

0.211 

e 


0.265 

0.031 


in which we are interested here, are Jupiter-mass planets embed¬ 
ded in a 2:1 MMR, while d and e are much less massive. A small 
amplitude of libration is observed (~ 5°, see jCorreia et al.|2()TQ ) 
for the 2:1 resonance between GJ 876b, c. The forced eccentric¬ 
ity ratio is 


^2/e// 


- ^6.5. 


(69) 


5.2. GJ 876b, c: 2:1 MMR, small amplitude of libration 


GJ 876 is a M-dwarf hosting four planets ( 

Delfosse et al.|1998i 

[Marcy et al.|1998||2001l|Rivera et al.|2010 

). The planet b and c. 


Since the system is observed with a small amplitude of libration, 
the stability constraint gives (Eq. ( [50| )) 


-I 

V2)eU 


:42. 


(70) 
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As for HD 60532b, c, the condition of convergent migration is 
fulfilled (see Eq. ([54|)). The surface density profile exponent is 
(see Eq. ( [48] )) 

ySs « -5.2. (71) 


We obtain again a negative value that correspond to an inverted 


profile. We have (see Eqs. ( |4^ , ( [58] ), and ( |59| )) 

^ = 6.4xl0■^ (72) 

Cl = 0.81, (73) 

C2 = 22. (74) 

Using these values and r » 42, we obtain (see Eq. ( |60| ) 

K » 80, (75) 

which corresponds to an aspect ratio of (see Eq. ( |49| )) 

-«0.11. (76) 

a 


As for HD 60532b, c we performed numerical simulations 
with different values of r (and adjusted values of K given by 
Eq. ([60])). The semi-major axes are initially 2 and 3.5 AU (period 
ratio of about 2.3), the eccentricities are 0.001 with anti-aligned 
periastrons and coplanar orbits. The planets are initially at peri- 
astrons (zero anomalies). The evolution of the semi-major axes, 
the period ratio, the eccentricities, the eccentricity ratio, and the 
angles are shown in Eig. The transition between decreasing 
and increasing amplitude of libration happens around r 20 
(see Eig. [^, while our analytical criterion gives a value of 42. 
Taking this refined value for r, the condition for reproducing the 
observed system with type I migration reads 


K 

— « 0.13. 




-3.5, 

58, 


(77) 

(78) 

(79) 


The density profile still needs to be inverted < 0) and the 
overall concl usions are the same. 

Note that [Lee & Pealej ( |2002| ) studied capture scenarios for 
this system using a slightly different model for the migration 
and damping and did not observe any evolution of the libration 
amplitude in their simulations. The authors used constant semi¬ 
major axis (Ta,i) and eccentricity (Tep damping timescales for 
each planet. In our study w e followe d the prescriptions of jPa- 
paloizou & Larwood] ( |2000| ) (see also [Goldreich & Schlichting 


2014| ) and considered constant angular momentum and ec¬ 


centricity (Tep damping timesc ales. We replaced these prescrip¬ 
tions with [Lee & Pe'ale ( |20021 ) prescriptions for the disk-planet 
interactions in our analytical model (following the same scheme 
as described in Sect. We found that the amplitud e of li- 
bration does not evolve in this case (in agreement with [Lee ^ 
Peale|2002] simulations). This difference between both prescrip¬ 


tions has_Jmportant consequences since in the case of [Lee &] 
Peale| ( [200^ prescriptions two initially resonant planets will stay 


locked in resonance forever while with the prescriptions we used 
the amplitude of libration can increase and the system can escape 
from resonance. The main difference between both prescriptions 
comes from the fact that with [Lee & Pealej ( [2002] ) prescriptions, 
the eccentricity damping does not affect the semi-major axes, 
while in our model the eccentricity damping terms contribute 
to the semi-major axes evolution (see Eq. ([^). Disk-planet in¬ 
teraction models suggest that the semi-major axes evolution are 
indeed influenced by the eccentricity damping effect of the disk 
(see [Goldreich & Schlichting||20141 ). We thus follow these pre¬ 
scriptions in our study. 


Tabl e 3, Orbital parame ters of HD45364b,c used in this study (taken 
from jCorreia et al.|2009| ). The stellar mass is 0.82Mo. 


Parameter 

[unit] 

c 

b 

m 

[Mj] 

0.1872 

0.6579 

P 

[day] 

226.93 

342.85 

a 

[AU] 

0.6813 

0.8972 

e 


0.1684 

0.0974 


5.3. HD 45364b, c: 3:2 MMR, large amplitude of libration 


The star HD 45364 hosts two planets ( [Correia et al.||2009 ) 
bedded in a 3:2 MMR. The forced eccentricity ratio is 


em- 


(-) ^2.5. (80) 

V2K11 

A large am plitude of libration is observed (~ 70°, see [Correia 
[et al.|2009] ), thus, the stability constraint gives (Eq. ( [5T] )) 

T^(^) ^6.3. (81) 

The condition of convergent migration is fulfilled (see Eq. ([54])). 
The surface density profile exponent is (see Eq. ( [48] )) 

Pi: < -1.6. (82) 


We obtain again a negative value that corresponds to an inverted 
profile. We have (see Eqs. ( [4^ , ( [5^ , and ( [5^ ) 

6 = 6.0 X 10■^ 


Cl = 
C2 = 


0.93, 

2.3. 


Using these values and r > 6.3, we obtain (see Eq. ( [60] )) 
K > 9.4, 

which corresponds to an aspect ratio of (see Eq. ( [49] )) 

- < 0 . 33 . 


(83) 

(84) 

(85) 

( 86 ) 

(87) 


We performed numerical simulations with different values 
of T and K (given by Eq. ([^). The semi-major axes are ini¬ 
tially 10 and 14 AU (period ratio of about 1.65), the eccentrici¬ 
ties are 0.001 with anti-aligned periastrons and coplanar orbits. 
The planets are initially at periastrons (zero anomalies). The evo¬ 
lution of the semi-major axes, the period ratio, the eccentricities, 
the eccentricity ratio, and the angles are shown in Eig. [^ Ac¬ 
cording to our simulations, the amplitude of libration increases 
for T < 10 (transition between 5-15, see Eig.[^ which is compa¬ 
rable with our analytical result (r < 6.3). 

It may seem surprising that the amplitude of libration does 
not increase much more rapidly for r = 2 than for r = 5 (see 
Eig.[^. Indeed, our study shows that the smaller r is, the more 
unstable the resonant configuration is (Eq. @). However, the 
evolution of the amplitude of libration does not only depend 
on the ratio r = Tg 1 / 7^2 but also on the absolute values of 
these damping timescales (see Eq. ([^). In our simulations, we 
fixed the migration timescale for the outer planet (r^, 2 ) and var¬ 
ied the three other timescales: 1^,1 = Te ,2 = 

Te l = TTm, 2 lK. The damping timescales are thus much longer 
for the simulation with r = 2 (K = 3.5) than for r = 5 (i^ = 8), 
in order to reproduce the same equilibrium eccentricities. This 
tends to slow down the increasing of the amplitude of libration 
and compensates the acceleration provided by decreasing r. 
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Fig. 4. Same as Fig. [^but for HD 45364b, c. We used r = +oo, 15, 5, 2 with K = 15, 12, 8, 3.5 respectively for the four shown simulations (four 
columns). The amplitude of libration decreases for the first two simulations (r = +oo, 15) and increases for the last two (r = 5, 2). The value given 
by our analytical criterion for the transition between decreasing and increasing amplitude is r ~ 6.3. 


6. Conclusion 


We obtained a simple analytical criterion for the stability of the 
resonant configuration between two planets during their migra¬ 
tion in a protoplanetary disk. We used the simplified integrable 
model of mean-motion resonances that we developed in |Delisl^ 
al.| ( |2014| ), and modeled the dissipative effect of the disk on 
the planets by four distinct timescales: 7^,2 (migration of 

both plan ets), and Tgj, Te ^2 (da mping of both eccentricities). As 
shown by [Lee & Peale| ( |2002| ), migrating planets that are cap¬ 
tured in resonance have their eccentricities excited by the mi¬ 
gration forces of the disk. The eccentricities reach equilibrium 
values between the migration and damping forces. However, this 
equilibrium can be unstable, in which case the amplitude of li¬ 
bration in the resonance increases until the system crosses the 
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separatrix and escapes from resonance ( [Goldreich & Schlicht- 
ing|2014| ). We showed here that the equilibrium is stable on the 
condition that Te,ilTe ,2 > (^ 1 /^ 2 )^// (ratio of equilibrium eccen¬ 
tricities). For observed resonant systems, it is probable that the 
equilibrium was stable during the migration phase. Otherwise, 
the planets would have escaped from resonance. This result al¬ 
lows to put constraints on the damping forces undergone by the 
planets. For instance, using prescriptions for type I migration, 
we show that a locally inverted profile is needed for resonant 
systems for which the inner planet angular momentum deficit 
(AMD) is larger than the outer planet AMD. We applied our an¬ 
alytical criterion to HD 60532b, c (3:1 MMR), GJ 876b, c (2:1 
MMR), and HD 45364b, c (3:2 MMR). We showed for all stud¬ 
ied systems that if the planets had undergone type I migration, an 
inverted density profile would be required for the resonant con- 
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figuration to be stable. All these planets are Jupiter-mass planets 
and are thus believed to open a gap in the disk and undergo type 
II migration. Our results confirm that classical type I migration 
cannot reproduce the observed systems. 

Our model is very general and is not restricted to type I mi¬ 
gration. Considering a scenario of type II migration which is 
much more realistic for the studied systems, our model still gives 
constraints on the migration process and especially on the eccen¬ 
tricity damping undergone by each planet. However, we could 
not find a simple analytical prescription for type II migration in 
the literature and thus could not derive constraints on the disk 
properties in this case. Having analytical prescriptions for type 
II migration would allow a more detailed analysis of these sys¬ 
tems. 

It would also be very interesting in the future to study small 
planets in resonance (with precise enough determination of or¬ 
bital parameters to be sure of the resonant motion). Indeed, for 
small planets, a type I migration scenario is more realistic. In 
this case, a local inversion of the density profile (as needed for 
the three systems of this study) would be more surprising. 
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Appendix A: Renormalization of coordinates 

The renormalized variables are constructed by dividing all ac¬ 
tions by the following constant of motion (seelDelisle et al. 12012 

|2^ -- 

r=^Ai+A2, (A.1) 

k\ 


where A/ = Pi is the circular angular momentum of the 

planet i. Noting with a hat the initial actions, the renormalized 
ones are defined by 

A/ = 

A/ 
r ’ 

(A.2) 

Gi = 

Gi 

r ’ 

(A.3) 

li = 

li 

r’ 

(A.4) 

D = 

D 

Y' 

(A.5) 

S = 

s 

f’ 

(A. 6 ) 


Expressions Q and (j^ are straightforwardly derived from these 
definitions. 

Note that the Hamiltonian and the time also have to be renor¬ 
malized (see [Delisle et al.||2012[ |2014| ) in order to preserve 
Hamiltonian properties. However, in this study, we consider dis¬ 
sipative forces that act on the system on very long timescales. 
As long as the conservative timescale remains short compared to 
the dissipation timescale, the long term evolution of the system 
is well described by the mean effect of the dissipation over the 
conservative timescale. Therefore, the rescaling of this conser¬ 
vative timescale will not influence the long term evolution of the 
system. 


Appendix B: Reducing to an integrabie probiem 

In the general case, the motion of two planets in a ^ 2-^1 reso¬ 
nance is described by two degrees of freedom, i.e. both resonant 
angles 

CTi = —A2 - Ai-mu (B.l) 

q q 

^2 , ki 

0-2 = — A 2 - Al - V72, (H.2) 

q q 

and both actions h, h- Let us note Xi the complex cartesian co¬ 
ordinates associated to these action-angle coordinates 

Xi = (B.3) 

The simplifying assumption introduced in [Delisle et al.| ( |2014| ) 
allows to reduce this generally non-integrable problem to a sin¬ 
gle degree of freedom one (thus integrabie). The only remaining 
resonant angle is 6 and the associated action is D. Noting u the 


related complex cartesian coordinate 


u = 

(B.4) 

the simplifying assumption reads 


0 = cos V 2 - sin , 

(B.5) 

u = cos 0 e“^^^’"'^vi -f sin 0 e“^^^’"''x 2 , 

(B. 6 ) 


where 0 , ci^eih o- 2 ,eii are constant angles defined such that the 
libration center is directed toward u (see [Delisle et al.||2014| ). 
Equation ( |B. 6 | ) shows how the simplified one degree of freedom 
model mixes both initial degrees of freedom, and especially, how 
the resonant angle 6 mixes both initial resonant angles cri, 0 - 2 . 
Equation ( p^ directly results from Eq. ( |B.5| ). 


Article number, page 11 of 


12 
























A&A proofs: manuscript no. DCL 


Appendix C: Evolution of the parameter S under 
dissipation 

In this section we show how to compute the evolution of the 
parameter S (Eq. ( p^ ) under a dissipation affecting the semi¬ 
major axes and eccentricities of the planets. The evolution of the 
renormalized circular angular momenta only depends on (cr/djj. 
These renormalized quantities are constructed such that (see Ap¬ 
pendix 1^ and [Delisle et al.|2Q14| ) 


Ai Pi ^ mi j- 

— = -—— — Sa, 

A2 Pi ypi m2 

and 

y^Ai -r A2 = 1 . 

ki 

We deduce 


AiIj = 
M\d = 


A 1 A 2 a 
2 a 
k 2 A1A2 a 
ki 2 a 


The evolution of e is then straightforward (see Eq. Q) 
q A 1 A 2 a 


^\d = M\d M\d = 


ki 2 


a 


(C.l) 

(C.2) 

(C.3) 

(C.4) 

(C.5) 


The evolution of the renormalized deficit of angular momentum 
li is given by 


ii 

- 2 


li 

d 



(C.6) 


/i 


= 2 ^ 

A 2 a 

h 

d 


P 2 a 


(C.7) 


ii 


= 2 ^ 


ki Ai 

a 

= 2 ^ 

A 2 - 1 d 

h 

d 

^2 

d 

ki 2 

a 

d ^2 

A 2 a 


We thus deduce 


D 

D 


= cos 0 — 

h 


. 2 h 
+ sm 0 — 

d ^2 


21 cos^ 0 — 


A 2 - sim 0 a 


d 

• 2 ^2 
-r sm^ 0 — 

d ^2 


Einally, we have 
^\d = i>\d - ^\d 

= 2|cos^0 ^ 


a 


• 2 ^2 

+ sin^ (b — 

. ^2 


D 


A 2 - sin^ ^ a 
a 

q A 1 A 2 a 
ki 2 a 


D 


(C.8) 


(C.9) 


(C.IO) 
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